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telecommunication  systems  is  proposed  for  the  first  time  for  nonlinear  Dynamic  Range 
Compression  image  restoration  of  blurred  signals  embedded  in  very  high  noise  environment. 
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the  A-law/p-law  outperforms  image  restoration  based  on  the  well-established  image  restoration 
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deconvolution  orders  have  been  analyzed  using  the  nonlinear  transfonn  method. 
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Introduction 


Dynamic  Range  Compression/Expansion  known  as  companding  (compressing-expanding)  is  a 
well-established  principle  for  recovering  the  signal  embedded  in  high  noise.  Dynamic  Range 
Compression  accompanied  with  conventional  image  restoration  might  be  useful  in  recovering 
images  embedded  in  a  high  noise  environment.  Dynamic  Range  Compression/Expansion 
nonlinearity,  when  applied  to  a  noisy  signal,  improves  the  signal  to  noise  ratio  in  areas  where  the 
signal  is  low  compared  to  the  noise  and  reduces  the  SNR  in  areas  where  the  signal  is  higher  than 
the  noise  level.  This  principle  has  been  used  for  improving  the  quality  of  acoustic  signals  in  the 
50 ’s  and  is  extensively  used  for  noise  reduction  in  tape  recording  which  is  limited  by  “tape  hiss”, 
which  is  high  frequency  random  noise.  Noise  reduction  systems  like  the  “Dolby”  and  “dbx”  help 
to  solve  this  problem  by  pre-emphasizing  (compression)  the  high  frequencies  before  recording 
onto  tape  in  order  to  make  them  higher  in  amplitude  than  the  tape  hiss  noise  with  which  they 
compete  and  then  upon  playback,  a  matched  de-emphasis  filter  (expansion)  is  employed.  These 
noise  reduction  systems  are  amplitude  sensitive  so  that  only  soft  high-frequency  sounds  are 
emphasized.  And  the  matched  de-emphasis  filter  scales  the  high  frequencies  to  their  right 
position  with  the  other  recorded  signals1. 

Dynamic  Range  Compression/Expansion  technique  has  also  been  used  in  telecommunication 
systems  using  a  nonlinear  element  for  simultaneously  compressing  the  data  (reducing  the 
signal’s  dynamic  range  for  uniform  quantization  process),  enhancing  the  SNR  and  at  the  receiver 
end  expanding  the  data  through  the  inverse  of  the  same  nonlinear  element  ’  . 

The  use  of  Dynamic  Range  Compression  techniques  and  power-law  nonlinear  transfer  functions 
are  not  new  for  image  processing  purposes.  Techniques  such  as  contrast  stretching4  which 
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compress  the  data  (by  a  nonlinear  transfer  function)  below  a  threshold  value  to  darken  the  image 
for  higher  contrast,  image  enhancement  using  gray-level  nonlinear  transformation  functions  (like 
logarithmic  transformation4  6  and  power-law  transformation4)  and  simultaneous  gray-level 
compression  and  contrast  enhancement  or  multiplicative  to  additive  noise  conversion  using 
homomorphic  nonlinear  filtering  has  been  examined4'6. 

In  this  paper,  Dynamic  Range  Compression  image  deconvolution  via  both  the  A-law  and  the  p- 
law  encoders  is  introduced  for  the  first  time.  In  image  deconvolution  of  a  noisy  and  blurred 
signal,  dynamic  range  compression  performs  the  following  dual  functions:  (1)  Enhances  the 
high  frequencies  of  the  degraded  signal  relative  to  the  low  frequencies,  reducing  the  blur;  (2) 
Enhances  the  signal-to-noise  ratio  in  particular  at  high  frequencies  even  where  the  signal-to- 
noise  ratio  is  very  low,  further  improving  the  finer  details  in  the  image.  Such  an  approach 
overcomes  even  the  optimal  Wiener  filter4  limitations.  These  limitations  are  described  in  the  next 
section.  Several  forms  of  distortion  functions  with  tested  input  images  are  studied  and  simulated. 
The  distortion  functions  such  as  motion  distortion,  atmospheric  turbulence  and  misfocusing 
aberration  is  examined.  A  human  face  (Lena)  is  utilized  as  an  input  image.  The  joint  spectra 
orders  for  this  deconvolution  technique  are  analyzed  using  the  nonlinear  transfonn  method  :  both 
A-law  and  the  p-law  nonlinear  functions  are  firstly  expanded  to  their  power-law  nonlinear 
saturation  functions,  then  the  contribution  of  each  order  has  been  calculated  and  finally  all  the 
orders  contributions  have  been  combined. 

Inverse  and  Wiener  Filters 

Inverse  filtering4  is  the  simplest  linear  approach  for  restoring  the  distorted  image.  One  can 
compute  the  estimated  restored  image  by  dividing  the  Fourier  transform  of  the  degraded  image 


F(u,v) 


G(u,v ) 
H(u,v ) 


F(u,v ) 


N(u,v ) 
H(u,v ) 
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(1) 


G(u,  v)  by  the  Fourier  transform  of  the  aberration  function  H(u,  v)  as  stated  in  the  following 
equation: 

Where  G(u,  v)  is  the  Fourier  transform  of  the  distorted  image  plus  noise, 

G(u,v)  =  F(u,v)H(u,v )  +  N(u,v )  (2) 

Equation  1  tells  us  that  even  if  the  distortion  function  is  known,  the  undegraded  image  can  not  be 

recovered  because  of  having  noise  N(u,  v)  in  the  image.  Because  the  noise  is  random  and  the 
Fourier  transform  of  a  random  function  is  random  and  unknown,  therefore  N  (u,  v )/H  (u,  v)  ratio 
could  easily  dominate  the  estimated  recovered  image  F(u,v)  and  create  a  noisy  undegraded 
image  and  if  the  noise  level  is  high  the  recovered  image  may  not  be  visible.  The  Inverse  filter, 
for  most  of  blur  functions,  is  a  highpass  filter  and  because  noise  spreads  over  a  broad-band  of 
frequencies  therefore  after  Inverse  filtering  the  noise  will  be  enhanced  in  high  frequencies.  To 
overcome  this  problem  the  Wiener  filtering  method  comes  into  the  consideration.  This  nonlinear 
approach  incorporates  both  the  degradation  function  and  the  statistical  characteristics  of  noise 
into  the  restoration  process.  This  method  considers  images  and  noise  as  random  processes.  The 


H(u,v )  ~  F(u,v )  ~ 

V 

1 

H(u,v) 2  F(u,v) 2  +  N(u,v ) 2 

A 

H(u,v ) 

W(u,v)  = 


Wiener  filter  expression  could  be  shown  as7 


(3) 


2 

where  F[(u  ,v)  is  the  distortion  function,  H*(u,  v)  is  the  complex  conjugate  of  H(u,  v),  \H(u,  v)\  , 

| F(u,  v)\  and  \N(u,  v)\  are  the  power  spectrum  of  the  distortion  function,  the  undegraded  image 
and  the  noise  respectively. 

Accordingly,  the  Wiener  Filter  consists  of  two  filters:  (a)  the  proportionality  filter  (the  first 
tenn),  (b)  the  inverse  filter  (the  second  term).  The  inverse  filter,  as  it  was  discussed  before,  is  the 
best  filter  for  image  restoration  in  a  non-noisy  environment.  The  proportionality  filter  represents 
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the  ratio  of  the  clean  blur  signal  power  spectrum  to  the  total  power  spectrum.  For  a  non-noisy 
environment,  the  proportionality  filter  has  100%  transmitivity  of  the  input  signal,  while  in  a  very 
high  noise  this  filter  transmits  nearly  0%  of  the  input  signal.  Therefore  the  process  behind  the 
Wiener  filter  consists  of  two  stages:  first  the  noise  filtering  via  the  spectrum  proportionality 
filter,  and  second  the  signal  deblurring  via  the  inverse  filter.  In  this  process  for  regions  when 
SNR  is  low,  the  signal  in  this  part  of  the  spectrum  will  be  lost.  Therefore,  the  Wiener  filter  is  not 
effective  in  harsh  operating  conditions  such  as  imaging  on  a  rainy  day  or  in  a  sandstorm.  When 
the  noise  is  much  higher  than  the  signal  across  the  full  spectrum,  both  the  noise  and  the  signal 
are  removed.  However,  our  proposed  dynamic  range  compression  deconvolution  is  a  new  signal 
recovery  technique  that  can  retrieve  the  signal  in  a  high  noise  environment  while  outperfonning 
both  inverse  and  Wiener  filters. 

Architecture 

The  architecture  of  image  compression  nonlinear  joint  Fourier  processor  via  A-law/p-law  is 
shown  in  Figure  1.  This  figure  shows  a  joint  Fourier  processor  setup  which  uses  a  single  input 
plane  for  both  the  reference  (r)  and  the  signal  (s)  images  and  jointly  Fourier  transforms  the  input 
via  a  lens.  The  joint  Fourier  spectrum  then  is  captured  by  a  CCD  camera  which  acts  as  a  square- 
law  receiver  and  is  responsible  for  mixing  the  spectra  to  produce  the  cross  products.  The  output 
from  the  CCD  camera  is  sent  to  the  A-law/p-law  digital  or  analog  receiver.  After  nonlinear  A- 
law  or  p-law  Dynamic  Range  compression  deconvolution  the  output  is  sent  to  a  spatial  light 
modulator  (SLM),  the  output  from  SLM  is  then  Fourier  transfonned  via  a  lens  to  produce  the 
processed  output.  As  shown  explicitly  in  the  theoretical  section,  the  CCD  camera  is  responsible 
for  combining  the  Fourier  transforms  of  the  impulse  response  and  the  distorted  image  to 
compensate  for  the  phase  distortion.  And  then  the  A-law/p-law  nonlinear  transformation  is 
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responsible  for  enhancing  high  frequencies,  and  hence  noise  reduction  and  signal  to  noise  ratio 
enhancement.  The  input-output  nonlinear  transfer  functions  of  the  A-law  and  the  p-law  receivers 


are  respectively  defined  as2 


1  +  ln  A\E\ 

RE)  =  — — — —  sgn(F) 
1  +  ln  A 


ln(i^|£p 
ln(l  +  /u) 


Where 


(1) 

(2) 


sgn  (E)  =  \ 


1 

-1 


for  E  >  0 
for  E  <  0 


(3) 


and  E  is  the  input  information.  In  a  joint  Fourier  processor  E  is  always  larger  than  0  and  hence 
sgn  (E)  is  always  1 .  The  parameters  A  and  p  controls  the  amount  of  compression.  In  the  standard 
systems  A=87.6  and  p=255. 

Figures  2(A)  and  2(B)  show  the  input-output  transfer  functions  of  A-law  and  p-law  encoders 
defined  in  equations  1  and  2  respectively.  Three  different  values  of  A  and  p  has  been  selected  in 
order  to  show  the  changes  in  nonlinearity  performance  of  the  encoders.  The  values  of  A=87.56 
and  p=  255  are  the  standard  values  used  in  telecommunication  systems. 

Figure  3  shows  the  proposed  architecture  of  A-law/p-law  joint  Fourier  processor  for  Adaptive 
Optics  applications.  The  reference  and  the  signal  information  correspond  to  a  point  source  (guide 

o 

star  or  multi  guide  stars  )  and  the  object  of  interest  are  passing  through  the  distortion  medium 
which  produce  h  +  f*h  joint  image  where  “h”  is  the  impulse  response  of  the  distortion  medium,  f 
is  the  object,  “f*h”  is  the  distorted  image  and  asterisk  is  the  convolution  symbol.  This 
information  is  sent  to  a  spatial  light  modulator  and  is  processed  via  the  joint  Fourier  processor 
shown  in  figure  1 . 
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Theory 


3 

In  this  section  we  analyze  the  deconvolution  orders  using  the  nonlinear  transform  method  .  In  the 
joint  Fourier  processor  shown  in  figure  3,  described  in  prior  section,  the  joint  Fourier  transform 
energy  captured  by  the  CCD  camera  is  given  by 

E('’*.'’?)  =  F\R(.vI,vr)  +  S<y„vf  (4) 

4/ 

where  R  is  the  Fourier  transform  of  the  reference  infonnation  (r),  S  is  the  Fourier  transform  of 
the  signal  information  (s)  and  vx  and  vy  are  the  spatial  frequency  coordinates. 

Assume  an  input-output  nonlinear  transfer  function  of  a  processor,  f(E),  with  Fourier  transform 
of  F(co)  where  E  is  the  energy  and  co  is  its  complement  frequency  variable.  This  function  can  be 
demonstrated  in  its  summation  series  format  as: 

f(E)  =  YJHk  (vx ,  vv )  cos[2ky0vy  +  k(</>R  -  </>s )]  (5) 

k= 0 

where  yo  is  the  separation  factor  between  the  reference  and  signal  information,  Or,  Os  are  the 
reference  and  signal  phases  and  Hk  (vx,  vy)  is  the  weighting  factors  which  is  given  by 

Hk(vx,vv)=jEik^  F(co) cxp[/co{R2  +  .S2 )}/, (icoRSjdco  (6) 

sk  is  the  Neumann  factor  (ck  =  1  for  k=0  and  c/(=2  for  k>  0). 

Therefore  theoretically  as  long  as  the  F(co),  Fourier  transform,  of  an  energy  function  is  known, 
substituting  the  F(co)  in  Eq.  6  would  give  the  correlation9  or  deconvolution  weighting  orders. 
According  to  Eq.  6,  in  order  to  calculate  the  deconvolution  weighting  orders  of  a  nonlinear 
transfer  function,  the  Fourier  transform  of  the  input-output  nonlinear  transfer  function  is  needed. 
In  this  section  firstly  the  Fourier  transform  of  the  input-output  nonlinear  transfer  function  of  both 
p-law  and  A-law  defined  in  equations  1  and  2  are  calculated. 
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In  order  to  calculate  the  Fourier  transform  of  the  above  functions,  it  is  more  convenient  to  first 
expand  both  of  these  nonlinear  transfer  functions  into  their  power-law  expansion  expressions  as 
follows: 


f{E)  =  f(Gp 


1 

2n-\ 


4y SNn(GpE ) 


(7) 


where  Gp  is  a  generic  parameter  that  corresponds  to  either  p  or  A  and  Gc  is  a  generic  constant 
defined  as 


G  = 


1 0  for  //  -  law 

1 1  for  A  -  law 


f(Gp)  is  a  generic  constant  which  is  given  by: 


AG  J  = 


f(ju)=—— — 7  f°r  the  //  -  law  encoder 
ln(l  +  ju) 

f(A)  =  — - -  for  the  A -law  encoder 

1  +  In  A 


Sm  is  the  generic  negative  nonlinear  saturation  function  defined  as: 
SngAGpE)  =  \-[spg(GpE)\ 


(8) 


(9) 


(10) 


Where  SpGp(GpE)  is  a  positive  nonlinear  saturation  function  which  is  defined  for  the  p-law  and 
A-law  encoders  respectively  as: 


SPGp  ( GpE )  = 


Spp(/JE)-- 

SPA  (AE) : 


liE 

2  +  fjE 
AE-l 
AE  + 1 


for  the  /j  -  law  encoder 
for  the  A-law  encoder 


(ID 


Each  order  of  the  negative  saturation  function  can  be  fractionally  expanded  to  its  first  order 
expansion  as: 


J  NGp  (2»+l) 


y  £■—»() 


Gnn 


1  AE  +  y  +  ns 


(12) 
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,,,,  2  for  the  ju-Iaw  encoder 

Where  y  =  \ 

1  for  the  A  -  law  encoder 


By  combining  all  equations  and  knowing  the  Fourier  transfonn  of  a  negative  function  as: 


F(co)  =  -i  sgn(/y)-^-cxp  i  c\p(icoEt) ) 


And  substituting  F(co)  in  equation  3,  the  weighting  orders  for  the  first  and  third  order  of  the  A- 
law  power-law  expansion  respectively  are: 


HklVx’Vy)  = 


{(1  +  AJa_  )(1  +  AJa+  T  y  1  +  AJI+  +  {(1  +  AJa_  )(1  +  AJa+  )}1/2 


Hk3[Vx’Vy)=HkM’Vy)X 

,  /  JA-  ,  JA+  1  ,  3 ,  |1  +  AJa+  ,  ll  +  AJa 


1  +  AJa_  1  +  AJaJ  Hn{vx,v 


1 1  +  AJ A  p  +  AJa+ 


■A2JA+JAH20l{vx,vy)  (16) 


kH2{vx,vy)( 


HuWxWv 


1  +  AJa_  1  +  aja+ 


1 +  AJA+  1  1  +  AJa_  1 

11  +  ^-  P  +  AJa+  H01(vx,v 


And  the  weighting  factors  for  the  first  and  third  order  of  the  p-law  power-law  expansion 
respectively  are: 


H  (v  v  )  =  - _ ^ (17) 

{(2  +  fjJA_  fl  +  juJA+  )f2  [ 2  +  fdI+  +  fl  + f.JA_\2  + nf A+)j 

HJVx’Vy)=HJVx’Vy)X 

r  /  J A-  JA+  \  6kHm{yx,vy)(  1 2  +fjJA_  )2  +  fjJA+  2  2  i  \ 

AP  +  PA  (2  +  frf.n ) J  fiZv.)  (  P  +  Pj  P+P,-)  *  (18) 

J^_+^W2+  IE Zf  +  IE Zfl  ■  1 


H\\Vx>Vy 


2  +  2  +  //A 


1 2  +  a-  \1  +  fFA+  H0\vx,vy 


Where 


JA+=(R  +  S)2 
J,x  =R2  +S2 
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Therefore  having  known  each  order’s  weighting  factors  and  combining  all  together,  we  can 
calculate  the  total  weighting  factors  for  both  A-law  and  p-law  encoders  as: 


Hk{vx,vy)=f(Gp\( 


2I  — 

2  n  - 


— 4Y 

i  h 


2n  - 1 


2,-d  k.vj 


where  Gc=l  for  A-law  and  Gc=0  for  p-law.  Substituting  equation  20  in  equation  5  would  show 
the  conversion  of  the  A-law  and  p-law  input-output  transfer  functions  defined  in  equations  1  and 
2  into  a  summation  series  as  follow: 


/(£)  =  £/(G, 


'  2n  - 1 


^  2/i- r 


2,-1  )ivX’vy)  cos [lky0vy  +  k(<f>R  -&)] 


For  more  details  on  calculating  the  above  equations,  see  the  Appendix. 

In  order  to  generate  the  gray  level  recovered  image,  the  above  equation  should  be  multiplied  to  a 
statistical  average  (generic)  filter,  Fs(vx,  vy)7.  The  A-law  and  p-law  joint  Fourier  processor  can  be 
used  either  for  correlation  and  convolution,  or  other  signal  processing  functionalities.  Before 
going  further  in  the  simulation,  the  order  expansion  in  equation  21  used  for  correlation  and 
convolution  is  discussed.  In  the  correlation  joint  Fourier  processor,  the  signal  is  expected  to 
include  the  reference  matching  template,  other  non-matching  signals,  clutter  and  noise  .  For  the 
most  simple  case,  where  the  reference  and  signal  information  are  the  same,  <Pr=  <&s,f  (E)  simply 
consists  of  a  sinusoidal  grating  with  periods  of  2kyo,  where  k  is  the  grating  period.  The  Fourier 
transform  of  equation  21  consists  of  a  periodic  delta  function  convolved  with  the  orders 
expansion  producing  the  appropriate  results  for  the  correlation  processing  functionality,  while  for 
the  convolution  processing  discussed  below,  only  the  first  order  expansion  is  suitable  to  perform 
the  functionality.  Assume  that  the  signal  Fourier  spectra  is  HF  and  the  distortion  Fourier 
transform  is  H.  In  the  deconvolution  before  the  compression  that  happens  in  CCD  camera  two 
terms  of  HFH*  and  H*F*H  are  produced.  These  terms  after  Dynamic  Range  Compression 
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produce  several  orders  in  equation  20.  The  signal  phase  Os=  Oh  +  Of  and  the  reference  phase  is 
0r=  &h-  For  image  restoration,  the  k  (Or-  Os)  =  ±kOp.  This  condition  is  satisfied  only  for  the 
first  order.  The  other  orders  have  phases  of  ±kOf,  which  produce  high  order  deconvolution  terms 
of  the  restored  image  convolved  with  coefficients.  The  higher  orders  can  clutter  the  restored 
images  therefore  it  is  important  to  know  the  energy  contribution  of  these  orders  compare  to  the 
first  order. 

Figure  4  shows  plots  for  the  first  three  order  joint  spectra  coefficient  (A)  for  A-law  (B)  for  the  p- 
law.  The  dashed  line  is  for  k=l,  the  dotted  line  is  for  k=2  and  the  dashed  dotted  line  is  for  k=3. 
For  small  A,  the  weighting  coefficient  value  is  higher  for  lower  orders.  This  is  evident  from  the 
plots,  at  A=10‘  the  first  order  coefficient  value  is  2.086  while  the  second  and  the  third  orders  are 
0.62.  The  reason  for  having  the  first  order  higher  than  the  others  at  small  values  of  A  is  that  the 
A-law  joint  Fourier  processor  is  almost  linear  for  small  A’s  (see  the  figure  3(A),  f(E)  as  a 
function  E).  Therefore  at  small  A  values  the  A-law  joint  Fourier  processor  functions  as  a 
conventional  processor  (described  in  chapter  1)  and  does  not  introduce  higher  orders. 

As  the  A  value  increases,  the  second  order  increases  to  higher  values  than  the  third  order.  All  the 
coefficient  are  maximized  at  the  same  value  of  A=  0.371.  After  the  maximum  point  up  to  almost 
A=3,  as  A  increases  the  coefficient  values  decrease.  After  the  plots  intersection,  the  third  order 
starts  increasing  to  the  higher  values  than  the  first  and  second  orders  while  the  second  order 
remains  the  intermediate  and  the  first  order  becomes  the  lowest. 

At  A=87.6,  the  first,  the  second  and  the  third  orders  values  are  1.199,  1.6  and  1.97  respectively. 
This  yields  that  the  signal  to  clutter  ratio  results  from  high  orders  interference  could  be  as  low  as 
0.2232.  Therefore  enough  separation  between  the  reference  and  the  signal  information  is 
essential  in  order  to  designate  separate  areas  for  the  orders  contributions. 
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The  p-law  coefficients  plots  illustrated  in  figure  4(B),  showed  both  similarity  and  dissimilarities 
to  the  A-law  coefficients.  The  p-law  plots  show  similar  results  as  the  A-law’s  plots  where  for 
small  A’s,  the  first  order  had  the  larger  value  than  the  second  and  the  third. 

Similar  to  the  A-law  plots,  the  second  order  coefficients  became  larger  than  the  third  order  in  the 
range  between  p=10'2  and  2.  In  contrast  to  the  A-law  plots  where  the  coefficients  take  their 
maximum  values,  the  p-law  coefficients  drop  continuously.  Similar  to  the  A-law  where  all 
coefficients  intersected  at  some  A  value,  the  p-law  coefficients  also  intersect  almost  at  p=5. 
Similar  to  the  A-law  plots,  after  the  intersection,  the  third  order  starts  increasing  to  the  higher 
values  than  the  first  and  second  orders,  the  second  order  remains  the  intermediate  and  the  first 
order  becomes  the  lowest. 

The  p-law  has  similar  behavior  as  the  A-law:  for  small  p  the  joint  Fourier  processor,  shown  in 
figure  3(B),  is  almost  linear  (f(E)  as  a  function  E)  and  hence  at  small  p  values  the  p-law  joint 
Fourier  processor  is  a  conventional  processor  and  does  not  introduce  higher  orders. 

Similar  to  the  A-law  at  the  standard  value  of  p  =255,  the  first,  the  second  and  the  third  orders 
values  are  2.268,  3.743  and  5.053  respectively.  This  yields  that  the  signal  to  clutter  ratio  results 
from  high  orders  interference  could  be  as  low  as  0.1301.  Therefore  for  the  p-law  also  enough 
separation  between  the  reference  and  the  signal  information  is  essential  in  order  to  designate 
separate  areas  for  the  orders  contributions. 

The  similarities  and  dissimilarities  in  the  coefficients  plots  of  the  p-law  and  A-law  are  attributed 
to  the  similarities  and  dissimilarities  in  the  negative  saturation  nonlinearity  of  the  A-law  and  the 
p-law  used  in  the  analysis.  Comparing  the  negative  saturation  nonlinearity  of  the  A-law  with  the 
negative  saturation  nonlinearity  of  the  p-law  shows  that  both  the  numerators  are  similar,  except 
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the  p-law  numerator  has  extra  nonlinear  term.  This  comparison  can  be  clearer  from  these 
equations  plots  shown  in  figures  5(A)  and  5(B)  for  A-law  and  p-law  respectively. 

Figure  5(A)  shows  a  secondary  maximum  for  the  A-law  negative  saturation  nonlinearity  around 
AE(vx,Vy)=1.5  which  could  be  the  reason  for  having  maximum  values  for  the  A-law  coefficients 
(figure  4(A)).  While  for  the  p-law  shown  in  figure  5(B)  the  negative  saturation  nonlinearity 
drops  continuously  and  hence  there  is  no  maximum  value  for  the  p-law  coefficients  as  it  was 
illustrated  in  figure  4(B). 

The  similarities  and  dissimilarities  in  the  coefficients  plots  of  the  A-law  and  p-law  from  their 
negative  saturation  nonlinearity  functions  point  of  view  were  discussed.  However  the  A-law 
coefficients  plots  for  A>1  values  (Fig.  4(A))  with  the  p-law  coefficients  plots  (Fig.  4(B)),  can  be 
considered  very  similar.  This  can  be  explained  from  the  A-law  and  p-law  input-output  nonlinear 
transfer  functions  plots  shown  in  figure  3  that  their  functionalities  are  almost  the  same  for  A>1 
and  p>0  values. 

Computer  Simulation 

Two  forms  of  computer  simulation  were  conducted:  direct  and  indirect.  In  the  direct  simulation 
both  of  the  noisy  blurred  information  and  the  aberration  impulse  response  were  set  in  a  zero 
array  with  a  separation  of  almost  256  pixels  from  their  centers.  Both  were  joint  Fourier 
transformed,  squared  and  nonlinearly  filtered.  In  the  indirect  deconvolution  simulation  only  the 
first  order  expansion  (k=l)  from  equation  21  was  used.  Other  orders  are  irrelevant  to  be  used  for 
deconvolution  purposes  however  for  correlation  application  these  orders  can  be  used.  Each  of 
these  simulations  represents  different  view  of  the  system  performance.  None  represent  what 
exactly  could  be  the  case  in  an  optical  system. 
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In  direct  simulation,  contributions  of  all  orders  can  lead  to  orders  aliasing;  while  only  the  first 
deconvolution  order  is  needed  to  simulate  the  deconvolution  through  an  optical  system.  Aliasing 
is  a  simulation  artifact  and  is  not  applicable  in  optical  systems;  in  particular  aliasing  from  high 
orders  which  interfere  with  the  output  result. 

In  indirect  simulation,  both  the  high  orders  aliasing  and  interference  contributions  are  removed. 
Since  the  aliasing  is  a  simulation  artifact,  it  can  be  removed  however  the  high  orders  interference 
exists  in  optical  systems.  Therefore,  none  of  these  forms  of  simulations  predicts  the  exact  optical 
system  results.  However  the  combination  of  direct  and  indirect  simulations  should  illustrate  a 
good  prediction  for  the  deconvolution  results. 

In  the  direct  simulation  a  joint  image  of  a  2  x  2  pixels  point  source  and  a  128x128  pixels  Lena 
face  were  located  in  a  1024  x  1024  zeros  (null)  array  while  they  are  separated  by  256  pixels 
center-to-center.  The  joint  input  is  convolved  with  an  aberration  function  to  generate  the  joint 
blurred  image  and  the  blur  impulse  response.  The  noise  is  created  by  a  random-number 
generator,  which  produces  white  Gaussian  noise  with  zero  mean  and  a  variance  of  1  and  is  added 
to  the  blurred  image.  Figure  6  shows  the  simulation  results  for  motion  aberration.  The  motion 
impulse  response  was  a  2  x  50  pixels  rectangle.  Figure  6(A)  is  the  joint  image  of  motion  impulse 
response  and  the  noisy  blurred  image  with  signal  to  noise  ratio  5,  6(A’)  is  the  gray  level 
recovered  image  using  the  A-law  deconvolution  for  A=87.6  and  6(A”)  is  the  gray  level 
recovered  image  using  the  p-law  deconvolution  for  p=255.  The  B  and  C  rows  are  the  same 
sequence  as  the  A  row  except  for  signal  to  noise  ratios  1  and  0.1  respectively. 

Figures  7,  8  and  9  are  all  the  same  sequence  as  the  figure  6  except  for  different  aberration 
functions  the  atmospheric  turbulence,  misfocusing  and  phase  respectively. 
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The  atmospheric  turbulence  was  defined  in  reference  10  as  ±l{y x,v v)  =  e  (22) 

where  P=0.0025.  The  misfocusing  impulse  response  was  considered  as  a  circle  with  a  diameter 
of  20  pixels  of  ones.  The  phase  aberration  impulse  response  was  a  uniformly  distributed  random 
phase-only  function  of  128  x  128  pixels. 

The  capability  of  the  proposed  system  in  image  recovery  without  high  orders  interference  was 
also  simulated  using  equations  15-18  substituted  in  equation  20.  In  the  direct  simulation  the  same 
128  x  128  pixels  Lena  face  is  located  in  a  512  x  512  zeros  (null)  array  as  an  input  image.  This 
input  is  blurred  using  the  same  blur  functions  described  in  indirect  simulation.  The  same  noise 
function  also  was  added  to  the  blurred  image. 

Figures  (10,  10’),  (11,  11’),  (12,  12’)  and  (13,  13’)  show  the  simulation  results  based  on  the  first 
order  nonlinear  transfonn  method  expansion  presented  in  equation  20. 

The  non-primed  figures  represent  the  results  of  deconvolution  for  the  following  blur  functions 
respectively:  motion,  atmospheric,  misfocusing  and  phase.  For  all  non-primed  figures  (A)  is  the 
noisy  blurred  input  with  SNR=  5,  (A’)  and  (A”)  are  the  corresponding  recovered  images  via  the 
A-law  and  the  p-law  decovolutions  respectively.  The  B  and  C  rows  are  sequentially  the  same  as 
the  A  row  for  SNR=1  and  SNR=0.1  respectively. 

The  primed  figures  represent  comparative  study  of  the  A-law  and  p-law  deconvolution  with  both 
of  the  Wiener  and  the  Inverse  filters.  In  all  of  these  figures  (A)  represents  the  recovered  image 
via  the  A-law  deconvolution  for  SNR=1  (the  (B’)  images  in  the  non-primed  figures),  (B)  is  the 
recovered  image  using  Wiener  filter,  with  the  expected  value  of  F,  (C)  is  the  recovered  image 
using  Inverse  filtering  and  (D)  is  the  recovered  image  using  Wiener  filter,  with  the  exact  value  of 
F.  As  it  is  evident  in  C’s  figures  (primed  figures),  the  Inverse  filter  was  incapable  of  either 
recovering  or  even  detecting  the  image.  In  these  results  noise  is  simply  covering  the  whole  array. 
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By  comparing  the  A’s  with  the  D’s,  it  is  clear  that  for  all  the  blur  functions  tested  here  the  image 
recovery  of  noisy  blurred  signal  with  SNR=1  via  the  A-law  deconvolution  was  always  better 
(with  a  generic  shape  of  the  low  pass  filter  or  the  expected  value  of  F)  than  the  Wiener  filter  with 
the  exact  value  of  F.  When  the  exact  value  of  F  was  replaced  by  a  generic  shape  of  the  low  pass 
filter  or  the  expected  value  of  F  like  what  it  was  used  in  the  A-law  results,  the  Wiener  filter  failed 
as  it  is  evident  in  B  for  all  the  primed  figures. 

Both  the  A-law  and  the  p-law  deconvolutions  results  showed  the  possibility  of  recovering  the 
noisy  blurred  image  for  all  common  forms  of  blur  functions  tested  here  with  SNR  as  low  as  0.1. 
They  are  associated  with  noise  conversion  from  low  frequencies  to  high  frequencies.  This 
frequency  conversion  leads  for  further  enhancement  of  the  signal  to  noise  ratio  through 
spreading  the  noise  over  the  entire  array.  While  in  image  recovery  via  the  Wiener  filter,  the  noise 
remains  nearly  concentrated  on  the  original  noisy  image  area. 

Conclusion 

This  paper  introduced  for  the  first  time  the  dynamic  range  compression  deconvolution.  Here  the 
dynamic  range  compression  deconvolution  implementation  via  both  the  A-law  and  the  p-law 
encoders  (used  in  telecommunication  systems)  was  proposed.  In  this  approach:  (a)  the  phase 
information  was  restored  via  the  combination  of  the  square  receiving  with  the  dynamic  range 
compression  and  (b)  the  amplitude  information  (gray  level  image)  was  restored  through  using 
lowpass  filtering. 

The  dynamic  range  compression  simultaneously  enhanced  the  signal  to  noise  ratio  as  well  as  the 
higher  frequencies  relative  to  the  lower  frequencies  which  leads  to  increase  noise  frequencies. 
The  deconvolution  theory  of  both  encoders  were  developed  using  the  nonlinear  transfer  method. 
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The  simulation  results  showed  better  performance  of  A-law  than  the  p-law.  Both  A-law  and  p- 
law  outperformed  the  Wiener  and  the  Inverse  filters. 

The  A-law  and  the  p-law  encoders  are  used  in  both  digital  and  analog  systems.  The  dynamic 
range  compression  severity  is  controlled  via  the  compression  parameters  A  and  p  for  A-law  and 
the  p-law  respectively.  The  simulation  results  were  improved  when  these  parameters  increased. 
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APPENDIX 


The  Input-Output  nonlinear  transfer  function  of  the  A-law  encoder  is 


g(x)  = 


1  +  In  (A  x  ) 
1  +  In  (A) 


where  in  joint  transform  processor,  x  corresponds  to 


E(vx  ,vv)=  R(  vx ,  vy )  +  S  (vx ,  vy ) 


R  is  the  Fourier  transform  of  the  reference  information  and  S  is  the  Fourier  transform  of  the 
signal  information.  Accordingly  equation  A 1  can  be  rewritten  as 


g(E)  = 


1  +  ln  (AE(vx,vy)) 
1  +  ln(^) 


The  nonlinear  transfer  function  of  A-law  encoder  can  be  expanded  into  summation  of  power-law 
saturation  functions 

g(E)  =  g{A )  1  +  2  (spA(AE(vx,vy))  +  UspA(AE(vx,vy)))3+UspA(AE(vx,vy)))5  +-1  (A4) 

V  J 


Where  g(A)  = - - - 

SV  7  1  +  ln  A 


AE(v  v  )-l 

(A4.1)  and  SPA(AE(vx,vy))  =  )  x  1  (A4.2) 


S  denoted  for  saturation,  P  denotes  for  positive  saturation  functionality  and  A  is  for  A-law. 
Alternatively  the  above  expansions  via  positive  saturation  nonlinearity  can  be  replaced  by 
expansions  that  consist  of  negative  saturation  nonlinearity  S  na  (AE): 


SNAn(AE(Vx,Vy)) 

2  n  - 1 


rf*-*  1+fe-25 


where  the  nth  order  of  negative  saturation  nonlinearity  is  defined  as: 

=  1  -  [spME(v„vr))] 
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The  first  three  terms  of  the  negative  saturation  nonlinearities  are 


$nai  (AE) 


2 

AE(vx,vy)  +  1 


f 

SnaME)  =  2 

V 


iAE(v„vy)f +\ 
(AE(V„vr)  +  lf 


A 


J 


SNas(AE)  =  2 


(vx ,  n,,  ))4  + 1  o(+4fi  (vx ,  vv  ))2  + 1 
{AE(vx,vy)  +  lJ 


(A7) 


(A8) 


(A9) 


In  the  most  general  form  each  of  the  negative  nonlinear  saturation  functions  can  be  expanded 
into  the  following  form: 


Wi)(^)=lii5 


a 


An 


AE  + 1  +  ns 


(A10) 


For  example  the  third  order  power-law  saturation  functions  it  can  be  expanded  in  the  following 
form: 


$NA3  (  AE) 


_  5\Ab ) 


(. AE  +  \f 


j.  sv  I  1 

(x  + 1  -  s\x  +  lX^  +  !  +  £■) 


£ — >0 


a 


a 


x  + 1  ■ 


-  +  - 


-  + - - — 

x  +  1  x  +  1  +  a 


(All) 


Where 


x  =  AE  (A  11.1) 

3  3  2 

—  (A11.2) 

lea" 


a 


A0 


(A  11.3) 


a 


4+i 


3  3  2 

- 1 - 1 - r 

2  as1 


(A  11.4) 
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Similarly  for  p-law,  the  most  general  form  each  of  the  negative  nonlinear  saturation  functions 
can  be  expanded  into  the  following  form: 


$ NM{2n+l)(AE)  ~  'Yj 


a 


/un 


AE  + 1  +  ns 


(A  12) 


For  example  the  third  order  power-law  saturation  functions  it  can  be  expanded  in  the  following 
form: 


SNui(AE)  = 


3(/jE)2  +6/uE  +  4 
(2  +  juEf 


3x  +  6x  +  4 


*'->0  (2  +  x  -  s\l  +  x\2 


=  lim 

+  x+s)  £->0 


- 1 - 1 - 

x+2-s  x+2  x+2+s 


(A13) 


Where 

x  =  juE 

(A  13.1) 

>  6  4 

au ,  -  3  —  +  — 

£  £ 

(A  13. 2) 

8 

(A13.3) 

-  6  4 

_  3  +  7  +  TT 

(A  13. 4) 

For  the  first  order  negative  saturation  input-output  transfer  function: 


RE)  = 


(A  14), 


l  +  Gp(E  +  E0 ) 

using  the  nonlinear  transfer  method3,  the  weighing  factor  for  /,-order  is  given  by 

HRk{V*’Vy)=~ 


-2  GpRS 


fi+Gp(/,_  +ea\i+gp(ja+  +i?0)]}1/2l  \+Gp(Ji++EA+{\+Gp{Ja_  +£0)][l+Gpf/,+  +£0)j}1/2 


(A15) 


Where  Gp  ;s  a  generic  parameter  which  could  be  either  the  beam  ratio  m  in  two-beam  coupling, 
the  A  or  //  in  A-Jaw  or  the  ju-law  power-law  expansions,  and  the  joint  spectra  momentum  are 
defined  as 
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(A  16) 


ja_  =  {r-s)2 
ja+={r  +  s)2 
JAx  =  RS 
JI+=R 2  +S2 


According  to  the  above  equation  the  zero  order  is 

Hr°  {[l  +  Gp  {J A_  +  E0  )Jl  +  Gp  (JA+  +  E0  )]}1/2 


(A  17) 


1 1  no  (vxvy)  corresponds  to  the  average  rectification  efficiency  of  the  joint  spectra  as  a  result  of 
processing  within  a  device  that  have  negative  saturation  nonlinearity.  And  also  the  k*  order 
power  which  corresponds  to  the  order  of  the  coefficient  amplitude  normalized  to  the  processor 
amplitude  rectification  efficiency  is  defined  as 


PKkk*’Vy)  = 


-2  GpRS 


1  +  Gp  ( JI+  +  £0)  + flj  +  Gp  (JA_  +  E0  )[l  +  Gp  (JA+  +  E  0 )]} 


.1/2 


(A  18) 


Then  equation  A15  can  be  represented  in  the  simplest  form  as 


H  Rk  [vx ,  v\  )  =  Hr  o  (vA. ,  vy  )PRk  (v_r ,  v  ) 


(A  19) 


For  the  small  incremental  change  in  the  value  of  A,  from  A  to  be  A/(l-s)  where  s  is  a  small 
number,  the  Hro  (vxyy)  the  rectification  efficiency  is  changed  to: 


H R0  (y x  ,  Vy  )  =  ' 


1 +—JA 
l±e 


1  +  —  J 
1  ±e 


~n  _1/2 


A+ 


[l  ±£1  A  +S~@A. _ 

J(1  +  AJa_)(1  +  AJa+) 


(A20) 


Where  I  a  and  O  a  defined  as: 


AJ 


A- 


AJ 


A+ 


'A  2(l  +  AJ  A_ )  2(l  +  AJa+  ) 


(A20.1) 


0  *JA+JA- 

A  4(1  +  AJaJi  +  AJa+) 


(A20.2) 
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Similarly,  the  coefficient  power  (the  k*  order  normalized  coefficients)  can  be  approximated  as: 


1  +  I  +  AJa_  |  1  ll  +  AJa+ 

L  .v  2  V  +  AJa+  2  P  +  AJa_ 
=  W[l  =FJAe-QAe2] 


PRk(vx,vv,±s^{PA)k 


1  + 


ks 


_ ks2 _ 

4Pa^(\  +  AJa_)(\  +  AJa+) 


(A21) 


Where  JA,  Qa  and  PA  are  defined  respectively  as: 

j,-*-  (l+Lp±PZ+Lf±PP) 

M  2  Vl  +  AJm  2  \1  +  AJa_) 


(A21.1) 


Qa  = 


k 

AP^(\  + AJ  A_\l  + AJ  M) 


(A21.2) 


PA  = 


2AJ 


Ax 


i  +  aj„  +  J)TaJZY^aJ7) 


(A21.3) 


Combining  the  contribution  of  both  the  envelope  and  normalized  convolution  orders  yields  that 
that  the  weighted  convolution  orders  as  results  of  incremental  change  in  A  to  be: 


HRk\VX’Vy’±£)  = 


J{1  +  AJaJ1  +  AJa+) 


(l  +  dA  +  £20^j 


-  2AJ . 


AJ i  +  +  y]{l  +  AJA+)(l  +  AJA_)' 


(l  +  sJa-s2Qa)  (A22) 


which  can  be  rewritten  in  the  following  form 

H  Rk  (a  ,  vy  ,±e)  =  HRk  (vx ,  vv  ,0Xl  ±  si  A  +s2Oa\\  +  zJ  A- s2  Q, ) 

=  HRk  (vx ,  vy  ,o)[l  +  {I  A- J A  )s  +  (0A  -  Qa  -  JAIA  )s2  ]  (A23) 

=  HRk(vx,vy,0Xl  +  a1£  +  b1a2) 


Where  HRk  (vx,  vy,  0)  is  the  value  of  HRk  (vx,  vy,  ±  s)  when  e=0,  a ;  =Ia-Ja  and  b /  =Oa-Qa~  J.a  I  a- 


It  was  shown  in  Eq.  All,  the  expansion  of  the  third  order  negative  saturation  nonlinearity  as  the 
following  format: 


SNA3(AE)  = 


3  (AE)2  + 1 
(AE  + 1)3 


=  liml 

£->0 


3  3 

- +  - 


3  3  2 

—  +  —  +  — 


■  +  ■ 


AE  + 1 -s  AE  + 1  AE  + 1  +  s 


(A24) 
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If  the  contribution  of  the  HRk  (vx,  vy,±  s)  from  Eq.  A23  for  all  the  terms  is  substituted  in  above 
fractional  expansion,  the  contribution  of  the  third  order  negative  coefficient  is  given  by: 

2  A 


HRk3  k .  Vy  ,±e)  =  H  A  k .  vy  >-4  |-3+4]-4^a-  k  ,yv,  o)+H  ,  k  A,,+^f|+3+4 

\A  £  £  J  £  £  £  J 

Substituting  the  ///«  (yv,  vy,  ±  e)  values  from  equation  A23  yields  that 


(A25) 


{^~~  +  “7  V»  k’A  41  -  a,s  +  b,s 2 )-  XaRk  {vx,vy  ,o)+  ^  +  k  k’U  ^X1  +  +  V: 2 )  (A26) 

=  ^k.^42  +  6(^-^)+4(d,-a  -7,/  J] 

Substituting  all  the  parameters  within  the  above  equation  yields  that 

Ha3{vx,vy)  = 


-2AJ . 


{(1  +  AJa_\  1  +  A/J}1/2  {l  +  AJI+  +  {(1  +  AJaJ  1  +  A/Jf 


f  A-  ,  A+  ] 

3i(l  +  AJIaA  {(l  +  A/,_Xl  +  ^A+)}1/2) 

fl  +  A/,t  :  jl  +  AJ A_ 

+  i  +  ^aJ 

\  +  AJa_  yi  + AJa+_ 

(A27) 


a2JaJa- 
{i+ajaJi  +  aja 


k[\  +  AJlf  +  {{\  +  AJA_Xl  +  AjJ}'l2)(  (  JA_  JA+ 

2AJa„  {(l  +  AJa_  Xl  +  AJAa  )}1/2 


2+  l1  +  AJ A+  +  \l  +  AJ A-  |{(l  +  A/.,  )(l  + A/.,+  )}1/2 
l  +  A/,_  1  +  A/.J  Vl  +  AJ.  Vl  +  AJ,  11  A+,> 


Using  the  relationship  in  equation  F2  the  above  equation  can  be  simplified  to 


HRki^*’Vy)=HRkl^X’Vy)X 


2  +  3  A 


J  A-  _  +  _  JA+ 


\\  +  Aj  A_  1  +  AJ  A+ 


3kH  (vx,vy) 

^nk’U) 


i  +  aja 

1  +  AJ , 


1  +  AJa 

I  +  AJ . 


■ A2JA+JA_H\0l(vx,vy ) 


kHR0l{VX’Vy) 


HRu[yx^y) 


JA-  .  +  .  U+ 


1  +  A| 

v  vi  +  ^u_  l  +  ^Ut 


|l  +  ^,+  , 

Ii+aja_) 

\ 

1 

1 +aja_  i 

1  +  ajaJ 

(A28) 


For  the  analysis  of  third  and  high  order  power-law  negative  saturation  functionality  of  the  p-law, 
the  same  steps  as  in  the  analysis  of  the  A-law  was  proceeded.  The  deconvolution  orders  can  be 
approximated  as 

/  .  _  A* 

(A29) 


-2/a/,. 


{(2  +  /t/4_)(2+/i/,+)}'^  Ip +  /i/, +  +  42 +  /i/,_)(2 +  <“/,+  ) 


5+3// 

r  ,  ++  1 

1  3*(2  +  /i//++>/(2  +  /1/J_X2  +  /i/ 

AJ?  + 

1(2  +  +J  (2  +  /A/,+  f 

*(2+ A 

+,x 

r 
r  ja- 

A/,+  A(2  +  +J(2  +  +j)f,  , 

L  (2  +  A/,42  +  /t/J+)  /AjA^(2  +  MfA_l2  +  fJA J  [  ' 

l  2  +  pj A- 

2  +  fJA_  |  2  +  /a/,. 

2 +  /«/,+  t|2  +  /t/,. 


2+juJa_  \2  +  /aJ 
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The  Eq.  A29  can  be  simplified  to  the  following  fonnat 


HRk^x’Vy)=HRkSy*’Vy)* 

Ja-  ,  -h±  \  6kHmfY*’vr) 

(2  +  nJA.)  (2  +  fjJ  A+ )  J  HRU[vx,vy) 

kHlSv^vy) 

HRu{vx^y) 


5  +  3// 


2  + 


2  +  1<J a-  ,  2  +  MfA+ 

2  +  ///,,+  V  2  +  AT*-  y 


+  fi2JA+JA-H2ROi(vx 


1  + A 


J , 


V 


2  +  AA_ 

2  +  /z/*_  2  +  fjJ A+  ]l  2  +  fjJA_  ]l  2  +  /jJ A+  j 


2  + 


2  +  /jJ  A 


HRoSy„vy) 


v/j 


,vy)  (A30) 
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Figures  Captions: 

Figure  1.  Joint  Fourier  processor  via  A-law/p-law 

Figure  2.  (A)  A-law  input-output  transfer  function  for  A=l,  5  and  87.  56,  (B)  p-law  input-output 
transfer  function  p=  0,  5  and  255 

Figure  3.  Adaptive  Optics  imaging  system  for  image  restoration  using  A-law/p-law 

Figure  4.  Joint  spectra  coefficient  for  (A)  A-law  and  (B)  p-law 

Figure  5.  Third  order  negative  saturation  nonlinearity  for  (A)  A-law  and  (B)  p-law 

Figure  6.  Motion  aberration  simulation  results,  (A)  the  joint  image  of  motion  impulse  response 

and  the  noisy  blurred  image  with  signal  to  noise  ratio  5,  (A’)  the  gray  level  recovered  image 

using  the  A-law  deconvolution  for  A=87.6  and  (A”)  the  gray  level  recovered  image  using  the  p- 

law  deconvolution  for  p=255.  The  B  and  C  rows  are  the  same  sequence  as  the  A  row  except  for 

signal  to  noise  ratios  1  and  0.1  respectively. 

Figure  7.  Atmospheric  turbulence  simulation  results  (A)  the  joint  image  of  motion  impulse 
response  and  the  noisy  blurred  image  with  signal  to  noise  ratio  5,  (A’)  the  gray  level  recovered 
image  using  the  A-law  deconvolution  for  A=87.6  and  (A”)  the  gray  level  recovered  image  using 
the  p-law  deconvolution  for  p=255.  The  B  and  C  rows  are  the  same  sequence  as  the  A  row 
except  for  signal  to  noise  ratios  1  and  0.1  respectively. 

Figure  8.  Misfocusing  aberration  simulation  results  (A)  the  joint  image  of  misfocusing  impulse 
response  and  the  noisy  blurred  image  with  signal  to  noise  ratio  5,  (A’)  the  gray  level  recovered 
image  using  the  A-law  deconvolution  for  A=87.6  and  (A”)  the  gray  level  recovered  image  using 
the  p-law  deconvolution  for  p=255.  The  B  and  C  rows  are  the  same  sequence  as  the  A  row 
except  for  signal  to  noise  ratios  1  and  0.1  respectively. 
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Figure  9.  Phase-only  aberration  simulation  results  (A)  the  joint  image  of  Phase-only  impulse 
response  and  the  noisy  blurred  image  with  signal  to  noise  ratio  5,  (A’)  the  gray  level  recovered 
image  using  the  A-law  deconvolution  for  A=87.6  and  (A”)  the  gray  level  recovered  image  using 
the  p-law  deconvolution  for  p=255.  The  B  and  C  rows  are  the  same  sequence  as  the  A  row 
except  for  signal  to  noise  ratios  1  and  0.1  respectively. 

Figure  10.  Motion  aberration  indirect  simulation  results  (A)  the  noisy  blurred  input  with  SNR= 
5,  (A’)  and  (A”)  are  the  corresponding  recovered  images  via  the  A-law  and  the  p-law 
decovolutions  respectively.  The  B  and  C  rows  are  sequentially  the  same  as  the  A  row  for  SNR=1 
and  SNR=0.1  respectively. 

Figure  10’.  A-law  deconvolution  result  for  the  motion  aberration  compare  to  Wiener  and 
Inverese  filters  result  (A)  the  recovered  image  via  the  A-law  deconvolution  for  SNR=1,  (B)  the 
recovered  image  using  Wiener  filter,  with  the  expected  value  of  F,  (C)  the  recovered  image  using 
Inverse  filtering  and  (D)  is  the  recovered  image  using  Wiener  filter,  with  the  exact  value  of  F 
Figure  11.  Atmospheric  turbulence  indirect  simulation  results  (A)  the  noisy  blurred  input  with 
SNR=  5,  (A’)  and  (A”)  are  the  corresponding  recovered  images  via  the  A-law  and  the  p-law 
decovolutions  respectively.  The  B  and  C  rows  are  sequentially  the  same  as  the  A  row  for  SNR=1 
and  SNR=0.1  respectively. 

Figure  11’.  A-law  deconvolution  result  for  the  Atmospheric  turbulence  compare  to  Wiener  and 
Inverse  filters  result  (A)  the  recovered  image  via  the  A-law  deconvolution  for  SNR=1,  (B)  the 
recovered  image  using  Wiener  filter,  with  the  expected  value  of  F,  (C)  the  recovered  image  using 
Inverse  filtering  and  (D)  is  the  recovered  image  using  Wiener  filter,  with  the  exact  value  of  F 
Figure  12.  Misfocusing  aberration  indirect  simulation  results  (A)  the  noisy  blurred  input  with 
SNR=  5,  (A’)  and  (A”)  are  the  corresponding  recovered  images  via  the  A-law  and  the  p-law 
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decovolutions  respectively.  The  B  and  C  rows  are  sequentially  the  same  as  the  A  row  for  SNR=1 
and  SNR=0.1  respectively. 

Figure  12’.  A-law  deconvolution  result  for  the  Misfocusing  aberration  compare  to  Wiener  and 
Inverese  filters  result  (A)  the  recovered  image  via  the  A-law  deconvolution  for  SNR=1,  (B)  the 
recovered  image  using  Wiener  filter,  with  the  expected  value  of  F,  (C)  the  recovered  image  using 
Inverse  filtering  and  (D)  is  the  recovered  image  using  Wiener  filter,  with  the  exact  value  of  F 
Figure  13.  Phase-only  aberration  indirect  simulation  results  (A)  the  noisy  blurred  input  with 
SNR=  5,  (A’)  and  (A”)  are  the  corresponding  recovered  images  via  the  A-law  and  the  u-law 
decovolutions  respectively.  The  B  and  C  rows  are  sequentially  the  same  as  the  A  row  for  SNR=1 
and  SNR=0.1  respectively. 

Figure  13’.  A-law  deconvolution  result  for  the  Phase-only  aberration  compare  to  Wiener  and 
Inverse  filters  result  (A)  the  recovered  image  via  the  A-law  deconvolution  for  SNR=1,  (B)  the 
recovered  image  using  Wiener  filter,  with  the  expected  value  of  F,  (C)  the  recovered  image  using 
Inverse  filtering  and  (D)  is  the  recovered  image  using  Wiener  filter,  with  the  exact  value  of  F 
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Figure  6 
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Figure  8 
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Figure  9 
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